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Abstract 

I introduce an innovative methodology for deriving numerical mod- 
els of systems of partial differential equations which exhibit the evo- 
lution of spatial patterns. The new approach directly produces a dis- 
cretisation for the evolution of the pattern amplitude, has the rigorous 
support of centre manifold theory at finite grid size h, and naturally 
incorporates physical boundaries. The results presented here for the 
Swift-Hohenberg equation suggest the approach will form a power- 
ful method in computationally exploring pattern selection in general. 
With the aid of computer algebra, the techniques may be applied to a 
wide variety of equations to derive numerical models that accurately 
and stably capture the dynamics including the influence of possibly 
forced boundaries. 
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1 Introduction 

The evolution of spatial patterns is an important area of research. One 
example of interest itself and serving as a model for other systems is Rayleigh- 



Benard convection ||19| , |T0|, |13|, e.g.] in which a fluid layer is heated from below 
and cooled from above. The flow flow self organises into an evolving pattern 
of upwelling and downwelling. One of the most useful models describing such 
pattern evolution, for example see ||, (9)], (2.11)] or [|lT], (13)], is the 
Ginzburg-Landau equation which we consider in in one spatial dimension: 

at = ra + cuxx — d\a\'^a (1) 

where a{x,t) is the (complex) amplitude of the pattern in any locale, sub- 
scripts denote partial derivatives, and r, c and d are specific constants. The 
Ginzburg-Landau equation describes the evolution of the complex amplitude 
of spatially periodic structures as they evolve and interact. The derivation 
of the Ginzburg-Landau equation in any specific pattern problem is funda- 
mentally based upon the underlying structure varying slowly in space-time. 
However, here in §|] we derive the discrete form 

CLj = raj + ^(«j+i — 2aj -|- aj_i) — d\aj\'^aj (2) 

of the Ginzburg-Landau equation ([l|) without ever invoking slow space-time 



variations. Following earlier research introducing holistic discretisation |117| , 
the derivation is rigorously based upon centre manifold theory 0, |, e.g.] and 
ensures the discretisation faithfully models the underlying system. 

One amazing consequence of this approach is that we straightforwardly 
incorporate physical boundaries into the discrete model. Boundaries have 
significant effects on the pattern evolution |T9| , ^ e.g.]. Previously there have 
only been limited, usually just linear p937-8, e.g.], arguments about how 
to incorporate such boundaries into the analysis to give boundary conditions 
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for the Ginzburg-Landau equation (|1]). But here we use nonhnear analyses 
to derive appropriately modified modifications of (|^) near the boundary to 
account for its effects. In we discuss time dependent Dirichlet and Neu- 
mann boundaries as two examples. Our approach resolves the subgrid fields 
within each element. Near a boundary the method resolves the modifications 
to the subgrid fields due to the infiuence of the boundary and thus naturally 
creates a discretisation appropriate to the specified boundary condition. 

As a simple prototype pattern evolution problem, we here consider the 
Swift-Hohenberg equation |2l|] in one spatial dimension: 



ut = ru-{l + d^yu-u-^ . (3) 

This equation has often been used, as we do here, to investigate issues in 
pattern selection [0, ||, e.g.], especially the infiuence of physical boundaries 
in 2D [jlO|, |l], e.g.]. Linearly, the Swift-Hohenberg equation (|^) has spatially 
periodic solutions u cx exp{Xt + ikx) for wavenumber k where the growth-rate 

X = r-{l-ey. (4) 

Thus for parameter r < all spatial modes decay, whereas for r > a band 
of modes near = 1 may grow. Exactly at the critical parameter r = 0, 
all modes decay except for the two neutral modes exp{±ix). We construct a 
holistic numerical discretisation based upon these neutral modes. 



2 Develop a centre manifold discretisation 

Here we develop a model of the system using invariant manifold theory. The 
analysis is essentially local in space, it resolves structures over a number of 
"convective rolls" . Previous use of invariant manifold theory ||^ has generally 
sought global models evolving on the so-called inertial manifold P2|, e.g.]. 



Exactly at criticality, long lasting solutions are 27r-periodic in x. As 
shown schematically in Figure |1|, introduce artificial internal boundaries every 
p periods, spacing h = p27i. These divide the domain into elements, each 
centred on a grid point Xj. Denote the field in the jth element by Uj{x,t). 
The internal boundary conditions (iBCs) to hold at each end of an element 
are chosen to be the nonlocal conditions 



duj 
^ dx 



(1-7) 



dx 



-h/2 
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Figure 1: schematics diagram showing a varying "roll" structure (solid curve) 
discretised by introducing artificial internal boundaries every period {p = 1) 
at odd multiples of vr (vertical bars). The "roll" field in the jth element is 
parametrised by the amplitude aj (discs). 



dx 




at X = Xj + h/2 . 



(5) 



- x=Xj-\~h/2 



8bti X 



h/2. 



(6) 



and the same for the second derivative v = Uxx ■ These iBCs are parametrised 
by 7 which controls the interaction and information fiow between adjacent 
elements. When this coupling parameter 7 = these iBCs reduce to condi- 
tions, 



u. 



dx 



-h/2 



Uj ± 



dx 



(7) 



. Xj+h/2 



requiring the solution in the j'th element is /i-periodic — the field and its 
derivatives at the right end of the element must match smoothly the field at 
the left end. Whereas when 7 = 1 these iBCs require continuity of the field 
and its derivatives between abutting ends of adjacent elements, 



dx 



Uj+l ± 



+1 



dx 



a.t X = Xj + h/2 . 



Thus when the coupling parameter 7 = each element is isolated from all 
others, whereas when 7 = 1 we ensure enough continuity between elements 
to recover the Swift- Hohenb erg equation (^ throughout the domain. We use 
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these iBCs to model the Swift-Hohenberg equation, when 7 = 1, by basing 
analysis on the discrete elements that are isolated at 7 = 0. 

Apply centre manifold theory based upon the linear dynamics in uncou- 
pled elements. By notionally adjoining the dynamically trivial equations 
7 = r = we treat all terms with a factor of the coupling parameter 7 
or the forcing parameter r as "nonlinear" perturbations, /i-periodic solu- 
tions of the linearised Swift-Hohenberg equation, ut = —(1 + then are 
exp(A„t + iknx) for any integer n with wavenumber /c„ = 2Tm/h = n/p and 
growth-rate A.„ = —(1 — k"^"^ . Thus within each element all modes decay 
exponentially except the two modes with wavenumber k±p = ±1 which are 
neutral. Hence "linearly," the solution in each element evolves exponentially 
quickly in time to the roll solution 

uj = aje'"' + bjC-'"' , (8) 

where aj and hj are the complex amplitudes of the rolls — the solution is real 
if and only if all hj are the complex conjugate of a^. Theory P, |], e.g.] then 
assures us that for the original nonlinear Swift-Hohenberg equation, coupled 
between elements by the iBCs d^-^), there exists a centre manifold Ai on 
which solutions are parametrised by the collection of evolving amplitudes a 
and h: 

M = Mj(a, 6, 7, r, x) (9) 
such that CLj = gj{a, b, 7, r) and bj = gj{a, b, 7, r) . (10) 

Moreover, theory assures us quite generally that solutions of the Swift-Ho- 
henberg equation exponentially quickly approach solutions of ([T0|). Lastly, 
theory asserts that the functions in (p|- |lO|) are determined by substitution 
into the Swift-Hohenberg equation and the internal boundary conditions and 
then solving to some order in the parameters. 

Computer algebra available from the author performs the tedious calcu- 
lations using an iterative algorithm []T3[. However, before undertaking the 



modelling we define the amplitudes precisely as the element averages 

X j-Xj+h/2 . I j-Xj+h/2 

(^iit) = T u{x,t)e~^'^dx and bJt) = - u{x,t)e^^'^dx . 

h Jxj-h/2 h Jxj-h/2 



Then the solution field is found to be 



6 aj — 2ifi6bj] + (4/i(5aj — 2i6 bj) x 



+ foj-e"*^ + ^e"*'' [6%j + 2ifi6aj^ + (AfiSbj + 2i6^aj^ x 
+ 0{^^ + + r^/^) , (12) 
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where A = \\a\\ + \\b\\ measures the overall amplitude of the solution fieldQ 
and where fi and 6 are central mean and difference operators, fiaj = (aj+i/2 + 
aj-i/2)/2 and 6aj = aj+1/2 — ai-1/2 respectively. These 0(7) terms show the 
leading order effect of neighbouring elements upon the subgrid field in each 
element; however, there is as yet no effect upon the evolution which is still 
neutral gj = (jj = C(7^). Using computer algebra to compute the next order 
of interactions and I find the evolution 

dj = raj + ^6^aj-3^^apj + Oi^^ + A^ + r^), (f3) 

bj = rb, + ^5%,-3^^ajb^^+Oi^^ + A^ + r^). (14) 

Evaluating these at the coupling parameter 7 = 1 we recover a numerical 
model for the Swift-Hohenberg equation which is just the discrete version (0) 
of the appropriate Ginzburg-Landau equation (p. 

In this approach there is little merit in discussing equivalent pdes to 
the derived discrete models (p!3|-[T^). The reason is that here the element 
size h must contain an integral number of rolls and so the limit h ^ is 
not physically valid. For the same reason, there appears to be no merit in 
seeking consistency with a pde as done for other holistic discretisations 



In this approach we should discuss the numerical model as it stands. But it 
is relevant to observe that the long- wave limit, slow variations in j, should 
and does reduce to the relevant Ginzburg-Landau equation (|1|). 



3 Boundary conditions are straightforwardly 
determined 

Now consider the boundaries to the physical domain. For simplicity suppose 
that the conditions at the boundary, say the left boundary, are either one of 
the two cases: 

u={-iyait) and m,, = (-If ; (15) 
or u, = {-iya{t) and u,,,= {-iy(3{t). (16) 

These physical boundary conditions are incorporated into the analysis by 
replacing the left-hand IBC (^) of the left-most element, say element j = 1, by 



a boundary condition corresponding to one of ([151) or ([Tq) . Consequently this 



A multinomial term 7 A™r" is 0(79 + A? + r«/2) if and only if / + m + 2n > q . 
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left physical boundary is, without loss of generality, located at x = Xi — h/2 . 
However, we want this left-most element, j = 1, to still have the neutral 
periodic solution (Q) when the inter-element coupling parameter 7 = 0. Thus 
we actually replace the left-hand IBC (@) of the left-most element by 



Ui - 



dui 
dx 



(1-7) 



Ui 



dx 



=F7 Ml 



dui 
dx 



x=xi+h/2 

± 2(-l)f7a(t) at X = xi - h/2, (17) 



and similarly for v = u^x- When the coupling parameter 7 = ([T7| ) reduces 
to requiring /i-periodic solutions. Whereas when 7 = 1 (|T^ reduces to 



:i±l)ni-(lTl) 



dui 
dx 



±2{-lfa atx = Xi~h/2, 



and similarly for v = Uxx, which for the upper choice of signs specifies the even 
derivative BC (|T5|) and for the lower choice of signs specifies the odd derivative 
BC (|16D. Other specific boundary conditions may be treated similarly, but 
here I just restrict attention to these two cases. 

Computer algebra constructs the centre manifold model for both alter- 
native boundary conditions: ([TED giving the upper alternative of plus/minus 
signs seen below; and ( [T^ ) giving the lower alternative. The leading order 
influence of the boundary forcing upon the field in the leftmost element is 



Ml = e'''ai + ^e'''[{-{2±i)ai + a2Tbi 
An 



2(±mi + 02 ± (1 ± 2i)hi - ih2)x\ 



, 7 a ix 
± -h— e 



± 



7^/3 



7 + 5z 

16 
3 + z 



3z 1 
— X H — 



+ e" 



16 
4/i 



I 1 - ; 

-X H 

4 96 



96 

{h^-l2x^ 



-{h'-12x' 



[(=Fai + ia2 - (2 =F i)bi + 62 



2(±(1 =F 2i)ai + ia2 =F ih + b2)x] 



± 



7-5i 

16 
3-i 



3i 1 
— X H — 



96 



'ih'-12x' 



16 



2 

-X + 



+ C(7^ + A^ + r 



3/2 



96 

a\ + \l3\) 



12x^ 



See there are two sorts of effects of the boundary: first, the subgrid field 
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Figure 2: spatial subgrid structure (|l^) of the coefficients of a (solid), j3 
(dashed) and its second derivative (dotted) when the boundary condition 
at X = Xi — IT is (|15D for an element size h = 2i[ . These give the field if all 
amplitudes = 6^ = . See the a and the second derivative curves are at 
the boundary x = about a value of one higher than the field in the bulk of 
the element — this is appropriate as a specifies the boundary value and (3 the 
second derivative. 

generated by the grid values of the leftmost two elements differs from that in 
an interior element (|12D; secondly, the forcing parametrised by a and (3 also 
adjusts the subgrid field. For example, in Figure ^ is plotted the basis of the 
field in the leftmost element when all grid values aj = bj = . See that they 
resolve a boundary layer structure on a scale of Ax ~ 1 — the resolution is 
crude because this is only the ffist approximation to the adjustment neces- 
sary to account for the boundary. By resolving the subgrid spatial structure 
near a boundary we will generate appropriate discretisations for the given 
boundary condition. The resolution of structure near the boundary is im- 
portant for deriving correct boundary conditions for models [^, p937, e.g.]. 
As demonstrated in Figure 0, our approach does some of the near boundary 
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analysis identified in as necessary for deriving boundary conditions for 
mathematical models such as the differential Ginzburg-Landau equation. 

The relevant evolution is obtained by analysing to quadratic terms in the 
coupling parameter 7. We find that 0,2 and 62 are identical to that obtained 
for the interior, ([T3|-p!^ with j = 2. However, the evolution for the amplitude 



of the rolls in the leftmost element adjacent to the boundary is 

472 2 7^ 

di = rai + — -(02 - 2ai =F 61) - 3ai6i =F — (1 - + /^) 

Al^ li 

+ 0(7' + A^ + rMa| + |/3|), (19) 

/I 2 2 

h = rh + ^{b2-2biTai)-3aiblT^{l + i){a + P) 

+ 0(7' + A^ + rM«| + |/3|). (20) 

Set the inter-element coupling parameter 7 = 1 to recover a discretisation. 
In the absence of boundary forcing, a = (3 = 0, we analyse a little of the 
dynamics near the boundary by supposing for illustrative purposes that ai = 
a2 = a and hi = b2 = a (for real solutions). The evolution equations (p!9|-|20|) 
then reduce to 

4 

d ~ ra — — (a ± a) — 3|apa . (21) 

First, consider the upper alternative when u and specified at the 

boundary, (plSj) . Linearly, the real part of a will decay at a rapid rate, a 
negative growth-rate of r — S/Zi^; whereas the imaginary part has a growth- 
rate r. Thus the amplitudes near the boundary will rapidly evolve to be pure 
imaginary. This corresponds to solutions u oc sin(x) as expected for the given 
boundary conditions. Secondly, consider the converse lower alternative when 
Ur and 

'^XXX 

specified at the boundary, (0). Linearly, the imaginary part 
of a will decay at a rapid rate, a negative growth-rate ofr — S/h"^; whereas the 
real part has a growth-rate r. Thus the amplitudes near the boundary will 
rapidly evolve to be purely real. This corresponds to solutions u oc cos(a;) 
as expected for the given boundary conditions. Thus the near boundary 
discretisation (p!9|-pOD selects for fields u with evenly spaced rolls located so 
that u is zero at the boundary in the case of specified even derivatives, and 
so that the derivatives of u are zero at the boundary in the case of specified 
odd derivatives. These properties agree with earlier more specific work |[T9| , 
e.g.] 

The presence of boundary forcing will push the system away from its 
otherwise natural equilibria. For example, for the upper alternative, (PT| ) 
breaks the symmetry in Q^(a) and predicts an equilibrium 3?(a) ~ —h{a+P)/8 
instead of zero and so there is a change in amplitude and phase of the rolls 
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near the boundary. Observe that time variations in the forcing, a and /?, are 
not significant at this level of approximation: this absence is surprising given 
the appearance of time derivatives of forcing in other boundary discretisations 



at finite grid size [|16|; the absence is linked to the definition (^) of the 
amplitudes because other amplitude definitions generate a and (3 terms in 
the model. Cox and Roberts 0] showed that special parametrisations could 
eliminate time dependence upon forcing in a centre manifold — evidently the 
amplitude definition (|1T]) at least approximates such a parametrisation. The 
important point is that our treatment of near boundary elements generates 
discretisations which naturally incorporate the forced boundary conditions. 



4 Conclusion 

We have introduced a rigorously based method for deriving discrete ampli- 
tude equations for the evolution of spatial patterns directly from the original 
PDEs. The method has been applied the the Swift-Hohenberg equation but 
the principles apply to a wide variety of PDEs governing spatio-temporal evo- 
lution. One remarkable advantage of this direct derivation is that we may 
also derive appropriate discretisations near any forced boundary. 

The generalisation to patterns in 2D space appears a straightforward 
generalisation of that employed for reaction-diffusion equations [|12[. Space 
would be tessellated in some regular manner by introducing appropriately 
periodic iBCs. Then a centre manifold model would be constructed with a 
finite number of basis wavevectors in each element analogous to the regular 
analysis of Skeldon et al . 

Computer algebra readily determines models of higher-order in the asymp- 
totic expansions, higher-order in both or either of nonlinearity or of stencil 
width. However, the IBCs used then here appear to break a symmetry 

of the Swift-Hohenberg equation. More research is needed into a better form 



of the IBCs as well as further applications of the methodology. 

Acknowledgement: this research is partially supported by a grant from 
the Australian Research Council. 



References 

[1] M. Bestehorn, R. Friedrich, and H. Haken. Travelling waves in nonequi- 
librium systems. Physica D, 37:295-299, 1989. 



AJ Roberts, February 1, 2008 



References 



11 



[2] J. Carr. Applications of centre manifold theory, volume 35 of Applied 
Math. Sci. Springer- Verlag, 1981. 

[3] J. Carr and R. G. Muncaster. The application of centre manifold theory 
to amplitude expansions. II. Infinite dimensional problems. J. DifJ. Eqns, 
50:280-288, 1983. 

[4] S. M. Cox and A. J. Roberts. Centre manifolds of forced dynamical 
systems. J. Austral. Math. Soc. 5, 32:401-436, 1991. 

[5] M. C. Cross. Boundary conditions on the envelope function of convective 
rolls close to onset. Phys Fluids, 25:936-941, 1982. 

[6] M. C. Cross, P. G. Daniels, P. C. Hohenberg, and E. D. Siggia. Phase 
winding solutions in a finite container above the convective threshold. 
J. Fluid Mech, 127:155-183, 1983. 

[7] M. C. Cross and A. C. Newell. Convection patterns in large aspect ratio 
systems. Phystca D, 10:299-328, 1984. 

[8] M. C. Cross, G. Tesauro, and H. S. Greenside. Wavenumber selection 
and persistent dynamics in models of convection. Physica D, 23:12-18, 
1986. 

[9] M. D. Graham, P. H. Steen, and E. S. Titi. Computational efficiency and 
approximate inertial manifolds for a Benard convection system. J. Non- 
linear Sci., 3:153-167, 1993. 

[10] H. S. Greenside and W. M. Coughran. Nonlinear pattern-formation near 
the onset of rayleigh-benard convection. Phys Rev A, 30:398-428, 1984. 

[11] L. Kramer and W. Pesch. Convection instabilities in nematic liquid 
crystals. Annu Rev Fluid Mech, 27:515-541, 1995. 

[12] T. MacKenzie and A. J. Roberts. The dynamics of reaction diffusion 
equations lead to a holistic discretisation. In R. L. May, G. F. Fitz- 
Gerald, and I. H. Grundy, editors, EM AC 2000 Proceedings. Proceed- 
ings of the fourth biennial Engineering Mathematics and Applications 
Conference, pages 199-202, 2000. 

[13] A. C. Newell, T. Passot, and J. Lega. Order parameter equations for 
patterns. Annu. Rev. Fluid Mech., 25:399-453, 1993. 

[14] A. J. Roberts. Boundary conditions for approximate differential equa- 
tions. J. Austral. Math. Soc. B, 34:54-80, 1992. 



AJ Roberts, February 1, 2008 



References 



12 



A. J. Roberts. Low- dimensional modelling of dynamics via computer 
algebra. Comput. Phys. Comm., 100:215-230, 1997. 

A. J. Roberts. Derive boundary conditions for holistic discretisations 



of burgers' equation. Technical report, fhttp : / /arXiv . org/ abs/math . 



NA/010622^ , 2001 



A. J. Roberts. Holistic discretisation ensures fidelity to Burgers' equa- 
tion. Applied Numerical Modelling, 37:371-396, 2001. 

A. J. Roberts. A holistic finite difference approach models linear dy- 
namics consistently. Mathematics of Computation, to appear, 2001. 

L. A. Segel. Distant side walls cause slow amplitude modulation of 
cellular convection. J. Fluid Mech, 38:203-224, 1969. 

A. C. Skeldon and M. Silber. New stability results for patterns in a 
model of long- wavelength convection. Physica D, 122:117-133, 1998. 

J. Swift and P. C. Hohenberg. Hydrodynamic fluctuations at the con- 
vective instability. Phys. Rev. A, 15:319-328, 1977. 

R. Temam. Inertial manifolds. Mathematical Intelligencer, 12:68-74, 
1990. 



AJ Roberts, February 1, 2008 



